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Abstract 

While from the energetic point of view supernova remnants are viable sources of Galactic cosmic rays (CRs), 
the issue of whether they can accelerate protons up to a few PeV remains unsolved. Here we discuss particle 
acceleration at the forward shock of supernovae, and discuss the possibility that the current of escaping particles 
may excite a non-resonant instability that in turn leads to the formation of resonant modes that confine particles 
close to the shock, thereby increasing the maximum energy. This mechanism is at work throughout the expansion 
of the supernova explosion, from the ejecta dominated (ED) phase to the Sedov-Taylor (ST) phase. The transition 
from one stage to the other reflects in a break in the spectrum of injected particles. Because of their higher 
explosion rate, we focus our work on type II SNe expanding in the slow, dense wind, produced by the red super¬ 
giant progenitor stars. When the explosion occurs in such winds, the transition between the ED and the ST 
phase is likely to take place within a few tens of years. The highest energies are reached at even earlier times, 
when, however, a small fraction of the mass of ejecta has been processed. As a result, the spectrum of accelerated 
particles shows a break in the slope, at an energy that is the maximum energy (Em) achieved at the beginning 
of the ST phase. Above this characteristic energy, the spectrum becomes steeper but remains a power law rather 
than developing an exponential cutoff. An exponential cut is eventually present at much higher energies but it 
does not have a phenomenological relevance. We show that for parameters typical of type II supernovae, Em for 
protons can easily reach values in the PeV range, confirming that type II SNRs are the best candidate sources for 
CRs at the knee. 

From the point of view of implications of this scenario on the measured particle spectra, we have tried to 
fit KASCADE-Grande, ARGO -YBJ and YACl-Tibet Array data with our model but we could not find any 
combination of the parameters that could explain all data sets. Indeed the recent measurement of the proton 
and helium spectra in the knee region, with the ARGO-YBJ and YACl-Tibet Array, has made the situation 
very confused. These measurements suggest that the knee in the light component is at ~ 650 TeV, appreciably 
below the knee in the overall spectrum. On one hand this finding would resolve the problem of reaching very high 
energies in supernovae, but on the other it would open a critical issue in the transition region between Galactic 
and extragalactic CRs. 
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1. Introduction 


Supernova remnants (SNRs) can provide the correct order of magnitude of energy injection in the form of 
accelerated particles to explain the density of CRs observed in the Galaxy (for reviews see, e.g., [18,8]). However, 
the issue of whether particles can be accelerated diffusively at the SN shock to sufficiently high energies is still 
open for debate. One thing that is clear is that effective magnetic field amplification is necessary, since in the 
absence of this phenomenon the maximum energy is bound to be in the ~ GeV region. This need was already 
recognized by [39,40], who found that the resonant instability induced by CR streaming could amplify the magnetic 
field to SBjB ~ 1 and allow the maximum energy to reach ~ 10^ — 10"^ GeV. This latter upper bound on the 
maximum energy, still well below the knee, derives from the limit 5B/B < 1, which is in turn related to the 
resonant nature of the instability considered [38,63,54,55,56]. Non resonant instabilities may in fact potentially 
lead to larger amplification factors, but on scales either much larger or much smaller than the particle gyroradius, 
so that larger fields would not necessarily reflect into increased particle scattering, which is what is needed to 
reduce the acceleration time and increase the maximum achievable energy. 

The pioneering papers of [41] and [12] showed how non-resonant modes can grow very fast ahead of a SNR shock. 
The fastest growing mode is characterized by wavenumber k much larger than the Larmor radius of the particles 
contributing the GR current, but the non-linear evolution of these modes is shown to form fluctuations on larger 
scales, up to the Larmor radius of the particles constituting the CR current. At that point, resonant scattering 
may in principle facilitate the return of particles to the shock surface and possibly increase the maximum energy. 
It has been argued that this non-resonant branch of the streaming instability may account for very high energy 
CRs in the Galaxy [13,52] and, in particular, may allow to reach the knee if applied to the case of supernovae 
expanding in the wind of their red super-giant (RSG) progenitor star. The role of winds of evolved pre-supernova 
stars for CR acceleration has been discussed several times in the literature, with different spins. Since massive stars 
often go through a Wolf-Rayet (WR) evolutionary phase and they are usually clustered within OB associations 
(superbubbles), several authors (see [15] for a recent review) have proposed that the collective effect of the WR fast 
rarefied winds may help accelerating CRs to very high energies and at the same time provide an explanation for 
some anomalies in the CR composition [33]. The most notable of such anomalies, and the one that is most solidly 
characterized from the observational point of view, is the ratio of abundances of two Ne isotopes, Ne 

[34], as observed in CRs, in comparison with its value in the solar environment. Additional constraints come from 
observations of in CRs, which impose limits on the rate of SN explosions in superbubbles, or, to state it 
more accurately, on the time when CR acceleration starts as compared with the time between two SN explosions 
in the same OB association. While this latter constraint is not a serious challenge for current models [15], it was 
recently pointed out ]44] that the average '^'^Ne/'^^Ne in superbubbles must approach the solar value, which would 
kill one of the main arguments in support of super bubbles as the main GCR sources. 

In this paper we focus on the explosion of SNe in the wind of the RSG progenitor star and we argue that these 
explosions are the best sites to reach PeV energies. A number of SNe expanding in their RSG wind have been 
detected, among which the best studied case is that of SN 1993J [51]. This source is especially interesting for our 
purposes because radio measurements provide us with direct estimates of the magnetic field strength at different 
times, through the detection of the synchrotron self-absorption feature [32]. This suggests a scenario in which not 
only the field is amplified at the SN shock, but the amplification turns out to be consistent with the expectations 
of the non-resonant instability [58]. Because of the high density of such winds, the highest energies are reached 
within a few decades after the SN explosion, namely before the beginning of the ST phase of the explosion. In 
this scenario, particles are accelerated from a pool that is mainly made of wind particles, rather than ejecta of 
other SN explosions. Eventually one may expect that the shock reaches the edge of the bubble excavated by the 
stellar wind, and at that point particle acceleration proceeds in the ordinary way, although the details will depend 
on whether the star was originally in an OB association or in the ISM. At such late times, the particles escaping 
the SN are less energetic and we will not be concerned much with this stage of the acceleration process. On other 
hand, it is worth recalling that most constraints on composition and composition anomalies are limited to such 
lower energies. 

The origin of the knee in the all-particle spectrum of GRs is inextricably connected with the issue of the end 
of the Galactic GR spectrum and the transition to extragalactic CRs: in Ref. [6], the authors try to model the 
transition region by summing the Galactic SNR contribution and the flux of nuclei of extragalactic origin, required 
to fit the Auger data [4]. The authors reach the conclusion that in order to satisfy the observational requirement 
of a predominantly light composition in the EeV region [4,2,3,60], an extra component of extragalactic protons 
is required. Such an extra component appears to be in good agreement with the proton spectrum as measured 
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by KASCADE-Grande [10]. However, both the proton and iron spectra measured by KASCADE-Grande in the 
energy region 10^® —10^® eV suggest that either the injection spectra are not cut off exponentially at the maximum 
energy (namely the number of particles decreases more slowly than an exponential at energies E ^ Em), or there 
is some, as yet unknown, class of sources with maximum energy much in excess of the knee energy. 

The observational situation has recently become even more confusing after the publication of the spectra of light 
GRs (protons and helium) by the ARGO-YBJ [26] and YAGl-Tibet Array ]36] experiments. Both measurements 
suggest a knee in the light component at ~ 650 TeV, appreciably below the knee. On one hand, it is clear that 
this finding would make the requirements in terms of particle acceleration less severe. On the other, however, 
it makes the requirements on the transition region much harder to accommodate. In other words, even in the 
scenario arising from ARGO data a population of sources able to reach maximum energy much above the knee is 
required to exist, thereby leaving the problem of the maximum energy little affected. 

In this paper we consider the viability of the two hypotheses above: we investigate whether non-exponential 
tails at E > Em are to be expected in the spectra of particles escaping supernova shocks, and we study the 
possibility of having supernovae that contribute maximum (proton) energies much higher than the knee. The 
paper is organised as follows: in § 2 we discuss the role of escaping particles for the estimate of the maximum 
energy achievable at a supernova shock. In § 3 we illustrate our calculations of the spectra of nuclei and in § 4 
we compare our predictions with the observed spectra, with particular emphasis on the knee and the transition 
region. In § 4.1 we discuss the implications of the recent measurements of the ARGO-YBJ and YAGl-Tibet Array, 
comparing them with KASGADE-Grande data. We conclude in § 5. 


2. Non-resonant instability excited by escaping CRs 


GR streaming leads to the excitation of both resonant Alfven waves and non-resonant, almost-purely-growing, 
modes ]12]. The former have a growth rate that, for typical parameters of a supernova remnant, leads to maximum 
energy <C PeV, mainly as a consequence of the saturation at SBjB ~ 1. On the other hand, the non-resonant 
instability is very fast growing, but the fastest growing mode has a wavelength much shorter than the Larmor radius 
of the particles generating the current, and also the wrong polarisation for resonant interaction with positively 
charged particles moving away from the shock surface ]7]. As a result, to first order, the non-resonant waves do 
not lead to efficient particle scattering, and would therefore be useless to the goal of increasing the maximum 
achievable energy. 

However, it was already suggested by ]12], and then confirmed by the numerical work of ]50] and more recently 
]22], that the non-linear evolution of these modes leads to the generation of power on larger spatial scales. Here 
it is useful to think in physical terms about the process that may be responsible for particle return to the shock 
from the upstream region. Let us consider a first generation of particles of energy equal to the current achievable 
maximum, E, moving from downstream to upstream. By definition, upstream there are no waves able to scatter 
resonantly such particles, hence they are going to move quasi-ballistically and escape the system at the speed of 
light (or fractions of it, depending on the level of anisotropy of the distribution of escaping particles). The current 
of particles at distance R from the explosion center as due to particles escaped when the shock location was at r 
can be written as: 

jcR{R, r) = ncR.r {R, EM{r)) eu,ft(r) = (1) 

with 

(EM/Eo)ln(EM/Eo) ^ = 0 
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1 + /3 
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\Em J 




( 2 ) 


In the above expressions we have used the relation between the number density of accelerated particles with 
energy E, ncR^E), at some location R upstream of the shock, and the energy density in accelerated particles 
at the shock ^CRp{''')v^{r), where ^cr is the GR acceleration efficiency. This relation is easily found through the 
transport equation for accelerated particles at a plane shock. The function 'L(Em) accounts for normalisation of 
the particle distribution function which is taken to be fs{E) oc and extending between a minimum energy 

Eo, that does not depend on time and a maximum energy Em which does depend on time, as we will see. In Eq. 
1, the factor {r/R)^ should account for the dilution due to spherical propagation of particles in ballistic motion. 
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The underlying assumption is that the large scale magnetic field is radial or absent. Whether this assumption is 
a realistic one in the two environments we consider in the following, namely the ISM or the wind of a massive 
progenitor star, may be questionable. In the latter case the zeroth order expectation is that the field should be 
toroidal, although one can envision that instabilities induced by the interaction between the slow and dense wind 
produced during the red supergiant phase of the progenitor and the faster and more dilute wind later blown as 
a Wolf Rayet star, could stretch the field lines and make them radial in at least a portion of the shock surface. 
In the case of a SNR expanding in the ISM, more caution should be taken, in that there is no reason why the 
upstream magnetic field should display a radial geometry. 

The fastest growing mode has a growth rate 


7 m = kMVA, 


( 3 ) 


where va is the Alfven speed in the unperturbed magnetic field Bq. The wavenumber where the growth is the 
fastest can also be easily estimated using the condition 


47r . 

Km Bo = — jcR, 
c 

which corresponds to balance between current and magnetic tension. The above expression implies that: 


fcMrn = (v)»l’ 


( 4 ) 


( 5 ) 


where rn is the particle gyroradius, ^cr is the CR efficiency and Va = the Alfven velocity. Upon saturation 
of the instability, the energy density in the form of amplified magnetic field in units of the ram pressure of the 
plasma can be estimated as (Vsh/ c){^cr/ Vi{EM/Eo)), that for typical values of the parameters is ~ 10”“* (for 
simplicity here we assumed that the spectrum of accelerated particles is ~ E~^). This reflects in about one order 
of magnitude larger pressure downstream of the shock (as due to compression), an estimate which is still a factor 
of 10 below that of [48], where the magnetic pressure downstream was assumed to be a fixed fraction (-^ 3.5%) of 
the ram pressure, independent of the shock speed. 

Numerical simulations show that the saturation of the instability occurs after Nt 5 e-folds [12], namely 

j 'yMdt = j(kMVA)dt « 5 . (6) 

we can rewrite this condition, at a point R upstream of the shock, 


Using the expression for kM derived in Eq. 4 
as [13,52,53,19]: 

47re ^CRp{r)Vsh{r)^ 


f 




dr ' 


(7) 

/o cEo ^^■kp{R)^{Em{R)) ' ’ 

where we have used the shock velocity to change the integration variable from time to shock coordinate. In Eq. 7, 
the integral in r is assumed to extend over all previous radii of the shock. In fact, since the values of the growth 
rate that are found for typical values of SN parameters are rather large (7 m Vsh/R), in case of ballistic motion 
of the escaping particles, one would expect that only a narrow range of radii close to R would give contribution 
to the particle current. Here we keep the form of Eq. 7 as given in Refs. [13,52,53], since the difference that this 
correction would introduce is only a factor of 2. 

In order to find an expression for the evolution of the maximum energy with time, we need to assume a density 
profile of the medium in which the supernova remnant is expanding. For a supernova expanding in the uniform 
ISM, one can assume that the ISM density is constant p{R) = pisM- For a type II supernova, instead, it is often 
the case that the explosion takes place in the wind produced by a red giant pre-supernova progenitor star. The 
density of the wind can be written as 


p{R) = 


M 


^-kRIv^ 


= p{Ro) ( ^ 


( 8 ) 


where M is the rate of mass loss of the red giant and Vw is the wind velocity. In the following we will simply 
write p{R) oc i?”™', having in mind the two cases of m = 0, corresponding to a uniform medium, and m = 2, 
appropriate for expansion in the progenitor wind. 

Differentiating Eq. 7 with respect to R, we obtain an implicit expression for the maximum energy (see Fig. 1): 

'&(Em(R)) = . (9) 
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3. From escape of accelerated particles to CRs 


3.1. Spectrum of the escaping particles 


For a supernova shock moving with velocity Vshit) in the surrounding medium (ISM or wind of the progenitor 
star) the shock radius is Rsh{t) = dt'vshit')- At each time t we assume that particles with given energy 
EM{t) = E can escape the accelerator. The number of particles that escape the shock at each give time is related 
to the current at the upstream radius R through: 


iVesc {E) dE = l^A-KR^dt 
e 


( 10 ) 


where A^esc (E) is the number of particles per unit energy, with energy E = EM{t), that are released in the ISM. 
Using the expression in Eq. 1 for the current one can rewrite: 


Ae.c (E) = 


f,CRPvl 


AtiR' 


dE 


( 11 ) 


with dR/dE = {dR/dA’){dA’/dE) to be computed using Eq. 2 and 9. In order to derive the first term on the right 
in the latter equation, we need to make explicit assumptions on the shock dynamics. For purposes of illustration 
we assume in the following a general power-law dependence of the shock radius on time: R oct^, which also implies 
Vs oc Using the latter dependencies in Eq. 9 it is straightforward to End for dR/d'^\ 

-1 
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R 


, ^-l 
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and rewrite Eq. 11 in a particularly useful way: 

Nssc{E) = 




-pvlR\{E) 


( 12 ) 


(13) 


EoS{X,m) 

where xi^) = {d/dE){l/'i’{E)) is a function of E alone, to be computed from Eq. 2. 

The advantage of Eq. 13 is that it makes it immediately clear that in the Sedov-Taylor phase, during which 
pVsR^ is constant with time, the CR spectrum released in the ISM is bound to be just given by xi^)^ hence 
only depend on the spectrum in the remnant, unless ^cr or Eq depend on time. It should also be noticed, however, 
that this dependence is not just a one to one identity and the spectrum of escaped CRs will not be the same as 
the spectrum of accelerated particles in the remnant, as already noted by [13],[22]. In fact, when computing xi^) 
explicitly from Eq. 2 we find: 
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(14) 


where the sign « refers to the fact that we have used the assumption E >> Eq. It is apparent that the power-law 
dependence of x{E) can only be —2 or steeper. Therefore the spectrum of CRs released during the Sedov-Taylor 
phase is the same as the spectrum of accelerated particles in the source only if the latter is E~^ or steeper, while 
for flatter source spectra NesciE) oc E~^. 

While in the Sedov-Taylor phase, only the environmental details are important while the ejecta density profile 
is irrelevant to the spectral slope, during the free expansion phase both density profiles are important in order 
to determine the position of the breaks and the maximum energy to which the released spectrum extends. This 
kind of information enters the calculation of the dependence of the shock radius on time, which we expressed as 
R oc t^. 

In the following, we take the value of A appropriate to describe the different evolutionary stages of the SNR 
from Ref. [24,45,59]. This work gives the remnant expansion law during the different stages for generic density 
profiles of both the SN ejecta pej oc R~^ and the ambient medium p oc R~^. The general expression for A is: 
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(15) 


where the ejecta profile appropriate to describe the two types of SN has k = 7 and k = 9 for type I and II 
explosions respectively. 

Once A is given, from Eq. 13 and 14 we can derive the spectrum released by the SN dnring the ED and ST 
phases as: 

iVe.c(E) oc pvjR^E-^^+‘> oc E -- (16) 

with e = 0 if fs{E) oc with P <9 and e = /? if /3 > 0. Inverting the time dependence of E as can be found 

from Eq. 9, and expressing the dynamical quantities as functions of E, one finally finds: 


NesciE) OC < 


Nesc[E) OC < 


E-(®+4') ED phase; 

m = 0, fc = 7 

£;-(2+£) 

E-(4+3") ED phase; 

m — 2,k = 9 


(17) 
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ST phase; 


This result is very interesting because it shows that the spectrum of escaping CRs integrated over time during 
the SN shock expansion is a broken power law, with an index close to 2 at low energies (belowEM) and steeper at 
higher energies. In more general terms, the spectrum of CRs released into the ISM is different from the postshock 
spectrum of accelerated particles, and is rather sensitive to how magnetic field amplification evolves in time [20,21]. 
A similar conclusion was previously obtained in the context of simple box models of CR acceleration at shocks 
[30,31]. At some energy ^ Em, the spectrum will eventually suffer an exponential cutoff. However such cutoff 
does not have an immediate phenomenological impact, because it would appear at energies at which the spectrum 
has already dropped appreciably due to the steeper power law. The transition between the two power laws at 
the maximum energy Em corresponds to the transition between the end of the ED phase, where most mass has 
already been processed, and the beginning of the ST (adiabatic) phase. Notice that in principle the highest value 
reached by particles at the shock is larger than Em at earlier times, but the flux of CRs with such energies is 
suppressed (though not exponentially but as a power law) because of the small amount of mass that is processed 
during the ED phase. 

In order to provide an estimate of the maximum achievable energy as a function of time, we proceed with 
modeling the transition from the ED to the ST phase of the SN expansion by using a parametric form for the 
time dependence of the shock radius: 


Rsh(t) — Ro 


to 


aXED 


- 1/a 


(19) 


where a is a smoothing parameter (in the following we use a = — 5 unless indicated otherwise), \ed and Xst are 
the indices describing the asymptotic time dependencies during the ED and ST phases, respectively, and Ro and 
to are the values of radius and time at the transition between the two phases. The value of Ro can be estimated 
by using the condition that at such distance the swept up mass equals that of the ejected material Mej : 


Mej = 47r / p{R)R dR. 

Jo 


( 20 ) 


where p{R) is the ISM density in the case of type I SNe and the density in the wind in the case of type II 
explosions. For the former type we then find, for the transition radius: 


Ro = 


3Mej 

47rp 


1/3 


f 7 nisM \ 3/3 
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( 21 ) 
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where p is the constant ISM density, which we take to correspond to 1 particle cm For type II events, instead: 

1(^11 , 1 ( -r 1 pc , (22) 
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The assumption of self-similarity for the time evolution of the SN shock leads to a generic form for the density 
evolution of the ejecta [59] 

Pej{R,t) = Ay. (23) 

where 


A = 


1 [3(fc - 3)Mejf 

4Trk [10(k - 5)EsJvf''^ 


10(fc — 5)Esn 


3{k - 3)Me, 


k/2 


(24) 


and Esn is the SN explosion energy. The condition that the density of the ejecta and the density of the wind be 
equal at the forward shock leads to the following expression for the time to'. 


to = 


D\ l/(k-m) 


(25) 


where B = pisM for type I and B = MjijpK'V.ui) for type II respectively. 
The shock velocity can be easily written as 
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(26) 


which can be used to calculate the maximum energy as a function of time, and hence the spectrum of escaping 
particles from an individual supernova. 

Before proceeding with detailed calculations, let us briefly discuss our choice of focusing on type II SNe as the 
best candidate producers of CRs at the knee. Let us consider the maximum energy that can be reached by each 
type of SN events, as given by Eq. 9 evaluated at the transition between the ED and the ST phase. To the goal 
of a simple estimate, let us specialise this equation to the case of a E~‘^ source spectrum. 

For type I events we can estimate: 


Em = 
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It is apparent that for standard parameters pertaining the two types of explosions, type II SNe can reach about 
one order of magnitude larger maximum energies than type la, and in particular that if particles are accelerated 
and escape as described here, they easily seem to be able to provide proton acceleration up to the knee. 

In order to illustrate how the maximum energy depends on the assumed spectral index within this framework, 
in Fig. 1 we plot the evolution with time of the instantaneous maximum energy a type II SNR can achieve for 
three different values of the source spectral index (all the other parameters have standard values, as specified in 
the caption). The figure shows that Em as derived from Eq. 9 and 2 changes by about 1 order of magnitude at 
the transition between ejecta dominated and Sedov phase (indicated by the vertical line in the plot) depending 
on the source spectral index: spectral indices steeper than 2 lead to lower values of Em compared to the estimate 
in Eq. 28, while the opposite is true for flatter spectra. 


3.2. CR Spectra observed at the Earth 

Here we adopt a simple model for the description of the transport of nuclei in the Galaxy. Assuming that 
the transport has reached a stationary regime and that it is dominated by diffusion, the equation describing the 
transport of a nuclear specie a is: 


Q I J-^a. Q j “h JaTT'0'sp,aV -— 0, 
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Figure 1. The maximum energy as a function of time after a type II SN explosion with EsN = 10^^ erg, ^CR — 0.1, Mej — Mq, 
M — MQyr~^ and Vw — lOfcm The different curves refer to different injection indices: p — —2 in red, p — —1.8 in green, 

p — —2.31. The vertical line identifies the time of beginning of the ST phase 


where Qa is the rate of injection per unit volume of CR of type a and asp,a is the cross section for spallation. In 
principle one should include an injection term of particles of type a due to the spallation of heavier nuclei, but 
as long as we focus on primary species and we are at sufficiently high energies, this contribution can be neglected 
[14]. Eq. 29 can be easily solved if we assume that both injection and spallation occur in a thin disc of size 2/i: 


n{z) = nd‘2h5{z) 




(30) 


Here the injection spectra Na{p) are calculated as discussed in Sec. 3.1, Ud is the gas density of the Galactic disc, 
assumed of thickness 2h and radius Rd- If we introduce the grammage 


X(p) = n4um,^, 


(31) 


and we impose a free escape boundary condition at \z\ = H, the solution of Eq. 29 in the Galactic disc can be 
easily written as follows: 

iV4p)» x(p) 

where p = 2hndmp is the surface density of the Galactic disc and Xa = mp/usp.a. Written in this form, the 
solution is very clear. At particle momenta for which the grammage is small compared with Xa the standard 
solution is recovered: 

iV„(p)K 


/o,c(p) = 


(33) 


2tvRIH D{p) ’ 

while in the range of momenta at which the solution is spallation-dominated, the observed spectrum reproduces 
the shape of the injection spectrum. In order to calculate the spectra of CRs at the Earth, an evaluation of the 
grammage is needed. In Ref. [47], the authors develop a leaky-box model in order to reproduce the GALPROP 
results and, for a source spectrum E~^ find a grammage in the form: 

- 0.6 


x{R) = igp-" 

where R is the rigidity of particles and f5 = v/c ' 
considering the reacceleration and find: 
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1 for 1 TeV particles. They calculate also the grammage 
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Given that the CR spectrum observed at Earth has an energy dependence as taken from TRACER and 

CREAM data [11,61], in our calculations we used a slope of the diffusion coefficient <5 = 2.65 — Pinj, where for 
Pinj we considered the values 2 and 2.31 as we will discuss in Sec 4. In the first case, the diffusion coefficient has a 
slightly stronger energy dependence than the above mentioned in the other case, instead, we use the exact 

form of Eq. 35. For the spallation cross section we adopt the simple formulation of Ref. [35]: 

asp{E)lmb] = (36) 

where A is the atomic mass number, and the two energy dependent coefficients are 

a{E) = 5.44 — 7.93 In 

/3(E) = 0.97 - 0.022 Zn 

It is worth noticing that one of the values of 5 we adopt here, and in fact the one that leads to most interesting 
results, is rather large, 5 = 0.65 (corresponding to pinj = 2), and this is known to lead to exceedingly large 
anisotropy in the CR signal at Earth [46,17]. 



4. Results 

Having shown that type II SNe are the most likely sources to accelerate protons up to the knee, in this section 
we focus on this type of explosions and attempt at fitting the all particle spectrum. To this aim we vary the 
different parameters, including the ejecta density profile. 

As discussed in § 3, the maximum energy up to which CRs can be accelerated in a SN is regulated by the escape 
of particles from the shock region and in the very early stages of the SN evolution, in principle, this energy can 
exceed the knee energy. However, the amount of mass processed at such times is relatively small and the net flux 
of particles is correspondingly small. The spectrum of accelerated particles that escape the shock region during 
the Sedov-Taylor phase and become CRs reflects the spectrum in the source if fs oc with /3 > 0, and is 

oc E~^ if /3 < 0. This spectral slope extends up to what we can name the effective maximum energy, Em, reached 
at the end of the ejecta dominated phase of expansion. We find that, at energies larger than Em, the spectrum is 
not falling exponentially, but rather as a power-law with a steeper index as compared to the lower energy trend. 
After escape from individual sources, CRs propagate diffusively through the Galaxy. At the high energies we 
are interested in here, the only relevant process during propagation is spallation, while other processes, such as 
reacceleration or Galactic winds, induce negligible effects and are therefore ignored in the present calculation. It 
is also worth stressing that the particles that are actually accelerated in the early phases of the SN explosion (say, 
within the first 100 years) are the highest energy particles; consequently, no direct constraint on these accelerators 
can be inferred from ^^Ni, which is often used to impose a lower limit on the time lag between the explosion 
and the acceleration phase. In any case, in this scenario, the material from which particles are accelerated to the 
highest energies is not polluted with SN ejecta, since it is basically the material expelled during the RSG phase 
of the pre-supernova star. 

As discussed in § 3, the value of Em is determined by the combination of the CR acceleration efficiency, ^cr, 
and the energetics of the SN, Esn- On the other hand, the flux of CRs at the Earth derives from a combination 
of CcH, Esn and the rate of SNe, K. The shaded (yellow) region in Fig. 2 shows the region of the ^cr — Esn 
plane for which Em > 1 PeV. The two panels refer to two different spectra at the source: E~^ is on the left and 
on the right. 

The need to fit the proton spectrum in the knee region imposes an additional constraint on the combination 
CchEsjvK. In Fig. 2 the lines indicate such a constraint projected on the ^cr — Esn plane for SN rates as 
specified in each panel. 

One can see that for the case of a ’’flat” source spectrum with /3 < 0, the standard SN energetics of 10®^ 

erg is compatible with the need to have maximum energy in the PeV range with totally reasonable values of the 
CR acceleration efficiency {^cr ^ 10%) and SN rate (K 1/30 yr“^). If the rate is decreased to 1/100 yr“^, higher 
energetics and higher efficiencies are required to reproduce the CR flux at Earth and still accelerate particles up 
to 1 PeV, whereas for a SN rate SR > 1/30yr~^, the knee simply cannot be reached: this is because for large SN 
rates the allowed CR energetics is low in order to avoid overproducing the flux at Earth. The conclusion is even 
more severe for SNRs with a steeper internal spectrum: in this case, the sources that can reach the knee and fit the 
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spectrum are extremely rare, with a frequency less than 1/800 yr“^; in addition the source energetics is required 
to exceed a few xlO®^ erg and the acceleration efficiency several tens of percent. 

The conclusion is that the most common type II SNe, expanding in their red supergiant wind, with standard 
energetics and an efficiency of order 10% are the most likely contributors of CRs up to the knee. In Fig. 3 we 




Figure 2. Left Panel: injection index p — —2. Right Panel: injection index p — —2.31. The shaded area in the ^CR ~ Esn plane 
indicates the allowed range of parameters to reach Em — 1 PeV. The lines indicate the combination of values of the parameters for 
which the observed proton flux is also fitted. Each line refers to a given SN rate (as labeled). 


show the spectra of the different nuclei, as well as the all-particle spectrum, for a population of type II SNe with 
a density profile of the ejecta described by A: = 9, an energy release corresponding to Esn = 10®^ erg, exploding 
at a rate = 1/30 yr“^. The spectra of nuclei are normalised to the data of CREAM for protons [61], PAMELA 
for Helium [5] and TRACER for CNO and Iron [11]. The diffusion coefficient is taken as described in Sec. 3.2. 
Clearly, in a simple approach such as the one discussed here, it is not possible to account for subtleties such as 
the harder spectrum of helium compared to hydrogen. We simply normalize the abundances in such a way that 
the CR composition at Earth is reproduced. We consider this assumption as a weak point of all models trying to 
explain the origin of CRs: there is still no comprehensive theory of acceleration of nuclei in SN shocks. In fact, it is 
known, and can be easily understood, that there is a preferential injection of heavier nuclei at a collisionless shock 
[27], but a quantitative theory of this phenomenon is not available yet. The situation is made even more complex 
by the fact that most refractory elements are embedded in dust grains, so that dust sputtering at shocks is bound 
to have en important role in determining chemical composition of CRs observed at the Earth. At present, the last 
and most comprehensive work on nuclei injection from dust sputtering is the one of Refs. [43,29], but a physical 
understanding, from first principles, of injection of particles with different masses and charges, is still lacking, so 
that one is forced to parametrize the injection of nuclei: this is what is commonly done and what we have done 
in this work, fixing the abundances so as to fit observations. 

Fig. 3 shows that for the values of the parameters that are typical of type II SNe the maximum energy of protons 
is in the PeV range, close to the knee at 3 x 10^® eV. The maximum energy of heavier nuclei is Z times larger. The 
spectrum of protons (as well as of other nuclei) is not cut off sharply at the highest achievable energy, but rather 
shows a severe steepening to a different power law. The slope of the high energy part is determined by the profile of 
the ejecta. For A: = 9 the spectrum at the source is steepened from E~^ to E~'^. The low energy part of the spectra 
of nuclei heavier than He shows a hardening that is due to spallation during propagation. The superposition of 
the spectra of all nuclei returns the all particle spectrum (black curve) that shows clear evidence for a knee, at 
energy of ~ 1 PeV. The high energy hardening of the protons spectrum is due to the introduction of an ad hoc 
extragalactic proton component with spectrum as found by the KASCADE-Grande experiment. 

In Fig. 4 we plot the maximum energy of accelerated particles as a function of time for an acceleration efficiency 
of 10% and different energetics of the parent SN, as labeled. One can see that more energetic SNe accelerate 
particles to higher energies and the onset of the Sedov phase in the wind occurs at earlier times. In any case, 
one should appreciate how in the framework of particle acceleration in SNe that explode in the wind of the pre¬ 
supernova star, the most important phase in the acceleration process occurs a few tens of years after the explosion, 
thereby implying a change of paradigm in which SNe, rather than SNRs, play a central role. The fact that in 
principle increasingly larger energies are reached at earlier times raises the issue of what is the minimum time 
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Figure 3. CR spectrum obtained with our model, fixing k — 9, Esn — 10^ erg and 9? = 1/30 yr~ . Protons(red), Helium (green), 
CNO (magenta). Iron (cyan), light component, p+He (blue) and the overall spectrum (black). At the highest energies we added 
by hand a constant extragalactic component. The proton knee is at E^ — 1 PeV, and the knee of the other elements is at 
E?,=Zx 


that we need to take into account in the calculation of the spectra of escaping particles. One obvious limitation 
conies from the fact that, after escape, particles need to cross the wind thereby suffering different kinds of energy 
losses. For protons, inelastic pp scattering occurs with a loss time scale Tpp — Ijnac, and it is unimportant when 


47 rr VuiTUp ^1^12 

-^- - S> - ^ r » 10 cm 


where reference values of the parameters have been adopted. For a shock velocity of km/s, this corresponds to 
times t S> 10® s. Given the large photon background in the early phases of the SN explosion, it is worth considering 



Figure 4. The maximum energy as a function of time after SN explosion with for different values of the total SN energy Esn- The 
change of slope identifies the transition from ejecta dominated phase to Sedov-Taylor phase. 


also inverse Compton scattering losses of protons, with rate of energy change: 
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Uph — 4 X 10 ty-p Eqpy GeV/ s 


where we approximate the photon energy density as Uph ~ Esm/{^R^). By requiring that E/{dE/dt) ^ r/c 
one gets that t ^ 6.3 x 10~^Ey^y yrs. For protons with PeV energy this constrains times useful for escape to 
later than the first ~ 2 x 10® s. For more energetic SNe (namely larger values of the shock velocity) this bound 
becomes less constraining. 

Photodisintegration of nuclei in the dense photon background may also become important. However, this is a 
threshold process, which is turned on when £ 7^1 > Eund ~ 10 MeV, Ebind being the binding energy of a nucleon 
inside the nucleus. This leads to a lower bound on the nucleus Lorentz factor 7 a > 10^. Above this threshold, the 
time scale for the process to become irrelevant, assuming a typical photon energy in the SN of e ~ 1 eV, is 

1 T 

- 3>->■ t S> 80 yrs. (38) 


The bound becomes less severe for more energetic SNe, namely for larger shock velocities. We stress once more 
that this bound is of some concern only for nuclei with rigidity quite above the rigidity of the knee. 


4.1. Comparison with high energy data 


Any understanding of the origin of CRs at the knee is bound to have important implications for higher energy 
CRs, and more specifically for the transition from Galactic to extragalactic CRs. In this section we discuss this 
issue in some detail. The standard lore about the origin of the knee in the all-particle spectrum is based on 
the premise that this feature coincides with a change of mass composition of CRs from light to heavy. In this 
perspective, the knee is associated with the end of the spectrum of the light component of CRs, say H and He. 
The natural implication of this line of thought is that the sources of Galactic CRs must be able to accelerate 
particles to at least the energy of the knee in the all-particle spectrum, say ~ 3 PeV. The important conclusion has 
stimulated endless discussion for the last thirty years on the sources and on the type of acceleration mechanism: 
in 1983, two seminal papers [39,40] reached the conclusion that even in the presence of resonant growth of Alfven 
waves the maximum energy reachable in SNRs could not exceed ~ 10® —10® GeV, thereby falling short of the knee 
by more than one order of magnitude. Measurement of the proton and He spectrum by KASCADE [9] seemed to 
confirm that the spectrum of the light CR component is steepening in the PeV region. A note of caution should 
be added, in that the information about mass composition is accessible to KASCADE only through a complex 
Monte Carlo dependent procedure. For instance using Sybill and QGSJet to simulate the showers, would lead to 
different spectra inferred for H and He. 

More recently, the KASCADE-Grande experiment [10] has measured the spectra and mass composition in the 
energy region E > 10^®'® eV: the spectrum of protons was claimed to show an-ankle like feature at energies 
E = iQi'i’ OSio os eV, where there is a hardening to a slope of about 2.7. This harder, high energy component was 
associated with the beginning of the extragalactic GR component. At the same time, the spectrum of Fe nuclei 
was observed to have a knee-\\\^e feature at E = eV, which also happens to be at ~ 26 times the 

energy of the knee. This feature was readily interpreted as the sign of a rigidity dependent knee in the spectra of 
the different chemicals. Notice however that the high energy part of the spectra (above the knee) does not appear 
to fall as an exponential, but rather as a steep power law. 

Qualitatively, this behavior seems similar to the one discussed above and connected with the transition from 
ejecta dominated to Sedov expansion of a SN shell. In Fig. 5 we show the spectra of different elements as calcu¬ 
lated using the model discussed above, and compared with the KASGADE-Grande data on protons and Fe. An 
extragalactic proton component with power law shape has been added to fit the highest energy data points 

of KASGADE-Grande. For the case k = 9, most appropriate for type II SNe (left panel in Fig. 5), a maximum 
energy Em — 3.7 PeV is required. Reaching this high energy requires a SN energy Esn = 2 x 10®^erg, a factor of a 
few larger than the standard Esn = 10 ®®erg, an efficiency ^cr = 20 % and a rate of explosion of SR = 1/110 yr“^. 
The value of Em is constrained by the need to fit the ankle-like feature on the proton spectrum as measured by 
KASCADE-Grande. For k = 7, less justified for this type of SNe, the high energy part of the injected spectra 
are somewhat harder and match the observations better. This is illustrated in the right panel of Fig. 5, where 
Em = 781 TeV and SR = 1/25 yr"®. 

Recent measurements of the spectrum of the light nuclei and the all-particle spectrum carried out with the 
ARGO experiment have raised serious doubts about the very foundations of the idea of light nuclei being accel¬ 
erated to ~ PeV energies [26,28]. ARGO finds a knee in the light component (H-fHe) at ~ 650 TeV, while the 
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Figure 5. Left) our best fit model for KASCADE-Grande data [10], with A: — 9, Esn — 2 X 10^^ erg, 0? — 1/110 yr~^, ~ 3.7 

PeV and ^CR — 20%. The highest energy data are fitted with an ”ad hoc” extragalactic component. Right) Our best fit model for 
KASCADE-Grande data [10], with A: = 7, Esn — 10^^ erg, 0? — 1/25 yr~^, E^ ~ 781 TeV and ^CR — 12%. 


all-particle spectrum is shown to have a knee at an energy compatible with all other measurements. This result 
appears to be compatible with some previous results obtained by experiments built at altitudes close to the shower 
maximum. The discrepancy between the ARGO and KASCADE results is an intrinsically observational issue and 
may be suggesting a rather poor knowledge of the development of showers and question our way to infer mass 
composition from shower-related observables. 

Since it is hard to imagine a way to make ARGO and KASCADE data compatible with each other, we decided 
to take ARGO data at face value and tried to infer the consequences that would follow. This is illustrated in Fig. 6, 
where we plot the results of our calculations and compare them with the ARGO data, shown as a shadowed area. 
The relatively large area derives from the fact that a few different analysis techniques have been adopted so as to 
establish the independence on the Monte Carlo procedure adopted of the fcnee-like feature in the light component. 

The best fit parameters we find in order to reproduce the ARGO data with our model are Esn = 10®^ erg, 
K = 1/15 yr~^ , Em = 507 TeV and ^cr '= 5.2%. While it is clear that the spectrum of light nuclei (H-|-He) is 
easily accounted for, it is equally clear that the all-particle spectrum knee cannot be reproduced within a simple 
approach such as the one discussed here. 
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Figure 6. Best fit model to ARGO data [26,28], with A — 9, Esn — 10^^ erg, 0f — 1/15 yr Em — 507 TeV and ^CR — 5.2%. 
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5. Summary 


Most of the models for the origin of Galactic CRs still rely on SNRs as being the most likely sources. Yet, the 
issue of whether SNRs are able to accelerate CR protons up to the energy of the knee, a few PeV, remains open. 
From the point of view of energetics, it seems that SNRs, in one flavour or another, are realistic candidates. From 
the point of view of the chemical composition, the anomalous ratio of ‘^'^Ne/'^^Ne observed in CRs has led to 
suggest that the bulk of Galactic CRs could be accelerated in OB associations [33,34,15], although this argument 
has been recently revised by [44] who reached the conclusion that OB associations cannot contribute the majority 
of CRs. 

In the present paper we focused on the issue of reaching the energy of the knee in particle acceleration during 
the expansion of a SN shock in the dense region occupied by the wind of a RSG, and stressing the physical aspects 
of the acceleration process, namely the nature of the instabilities that are required to bring the magnetic held 
strength and topology to the level that is suitable for particle acceleration. In particular, we investigated in detail 
the implications of the so called Non-Resonant Hybrid (NRH) instability described by [12,52,13] by computing 
the maximum energy and the overall particle spectrum produced during the whole SN expansion, both in the ED 
and ST phases. The general idea is that at any given time, the particles that escape the system with the highest 
available energy E produce the turbulence needed to scatter the next generation of particles with the same energy, 
thereby leading to acceleration to higher energies. We discussed the role of the NRH instability in both type la 
and type H SNe, in the context of three models of particle acceleration at the SN shock: 1) a flat spectrum oc E~^\ 
2) a hard spectrum oc with /? < 0, and 3) a steep spectrum oc f;“G+/ 3) /3 > 0. The maximum energy 

and the spectrum of escaping particles was calculated in each case. 

The effective maximum energy, which is reached at the beginning of the ST phase in all cases, is found to be 
at most of order a few hundred TeV for type la SNe exploding in the ISM (with the typical values of the relevant 
parameters). For type II SNe, the situation is quite more complex because these explosions often occur in the 
wind environment created by the pre-supernova red-giant star. Particle acceleration during the expansion of the 
SN shock in the dense slow wind of the RSG progenitor is found particularly promising from the point of view 
of exciting the NRH instability. In this case, the instantaneous maximum energy is shown to always decrease as 
a function of time after the explosion. However, the ST phase typically starts a few tens of years after the SN 
event and most of the mass is processed by the end of the ejecta dominated phase. These facts have important 
consequences. The early beginning of the ST phase allows for a maximum energy which is still high; the effective 
maximum energy does not coincide with a cut-off, but rather with a steepening of the spectrum, since higher 
energy particles are indeed produced during the ED phase. Our calculations show that a type II SN with standard 
energetics [Esn = 10®^ erg, ^cr = 10%) can accelerate CRs up to knee energies, Em ~ 1 PeV. 

It is interesting to appreciate how particle acceleration at PeV energies typically occurs ^ 10 — 30 years after 
the SN explosion, which represents a change of paradigm with respect to the standard SNR paradigm where the 
highest energies are reached several hundred years after explosion. In this scenario, the probability of catching a 
PeVatron in action in our own Galaxy (for instance through its gamma ray emission) is very low. 

For given parameters of the SN explosion, the maximum energy reached at the beginning of the ST phase is a 
function of the CR acceleration efficiency while the spectrum observed at the Earth also depends on the rate of 
SN explosions. We calculated the required efficiency and rates necessary to reach the knee and fit the overall CR 
flux: for steep spectra at the shock, the escape current is lower and this forces one to have larger SN energetics 
and larger CR acceleration efficiencies, which is counterintuitive for a steep spectrum of accelerated particles. For 
flat E~^ spectra at the shock, the spectrum of escaping particles is also flat and PeV energies can be reached 
for ordinary parameters of the SN explosion. Reaching PeV energies is even easier for hard spectra at the shock, 
however in both the flat and the hard case, one has to face the well known severe problem with CR anisotropy 
[46,17,49]. 

The comparison of the expected spectra of light nuclei with the recent data of KASCADE-Grande and ARGO 
appears very problematic: while individually both sets of data can be fitted with reasonable values of the SN 
parameters and CR acceleration efficiencies, there is no model that can fit both at the same time. KASCADE- 
Grande data require a SN energy Esn = 2 x 10®^ erg, an efficiency ^cr = 20% and a rate of explosion of SR = 
1/110 yr“^ that lead to a maximum energy Em — 3.7 PeV, whereas we can fit ARGO data with Esn = 10®^ erg, 
SR = 1/15 yr~^ and ^cr — 5.2% that lead to Em — 507 TeV. The disagreement between these two experiments 
is a purely experimental problem, that requires a serious and careful assessment of the systematic uncertainties 
involved in the adopted experimental techniques. 
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